# Libraries ----

library(dplyr)
library(tidyverse)
library(estimatr)
library(texreg)
library(ggplot2)
library(RColorBrewer)
library(sf)
library(tmap)
library(fect)
library(panelView)
library(patchwork)
library(DirectEffects)
library(mediation)

# Data ----

# load data for the outcome measure

df <- read.csv("democracy_measures_paper.csv")

# construct the treatment variable

df <- df %>%
  mutate(spending_ban = ifelse((st == "AK" | st == "AZ" | st == "MT" | st == "WY" | 
                                  st == "CO" | st == "TX" | st == "OK" | st == "SD" |
                                  st == "ND" | st == "MN" | st == "IA" | st == "WI" |
                                  st == "MI" | st == "OH" | st == "PA" | st == "WV" |
                                  st == "KY" | st == "TN" | st == "NC" | st == "NH" |
                                  st == "MA" | st == "RI" | st == "CT"), 1, 0))

df <- df %>%
  mutate(post_cu = ifelse((year >= 2010), 1, 0))

df <- df %>%
  mutate(treatment_dummy = ifelse((spending_ban == 1 & post_cu == 1), 1, 0))

## Figure 1: Treatment Status

us_shp <- st_read("cb_2018_us_state_20m/cb_2018_us_state_20m.shp")
us_shp <- us_shp %>%
  filter(STUSPS != "PR" & STUSPS != "DC") %>%
  left_join(df, by = c("NAME" = "state")) %>%
  tigris::shift_geometry()

map <- tm_shape(us_shp) + 
  tm_polygons(
    col = "spending_ban",
    palette = "Blues",
    n = 2,
    border.col = "white",
    title = "",
    labels = c("Control States", "Treated States")) +
  tm_layout(
    frame = FALSE)

map

# Descriptive Results ----

## Figure 2: Pre- and Post-Treatment Trends

plotdata <- df %>%
  group_by(spending_ban, year) %>%
  summarize(mean_democracy = mean(democracy_mcmc))

ggplot(data = plotdata, aes(x = year, y = mean_democracy,
                            group = spending_ban, color = factor(spending_ban))) +
  geom_line() +
  geom_vline(xintercept = 2010, linetype = "dotted") +
  scale_x_discrete(limits = c(2001, 2005, 2010, 2015, 2018)) +
  labs(title = "Citizens United and State Democratic Performance",
       x = "Year",
       y = "Average State Democracy Index",
       color = "Treatment") +
  scale_color_brewer(palette = "Paired", labels = c("No", "Yes")) +
  theme_minimal()

# Main Results ----

# main TWFE model

did_1 <- lm_robust(democracy_mcmc ~ treatment_dummy + factor(st) + factor(year), 
                   clusters = st, data = df, se_type = "stata")

## Table 1: Effect on State Democratic Performance (Bivariate Model)

texreg(list(did_1), include.ci = FALSE, 
       omit.coef = c('st|year|Intercept'), 
       custom.coef.names = c("ATT"),
       caption = "Effect on State Democratic Performance: Bivariate Model",
       custom.model.names = c("Model 1"),
       custom.gof.rows = list("State/Year Fixed Effects" = c("Yes")))

# FEct approach

set.seed(2010)

fect_result <- fect(democracy_mcmc ~ treatment_dummy, data = df, index = c("st", "year"),
                    method = "fe", force = "two-way", se = TRUE, vartype = "jackknife")

print(fect_result)

## Figure 3: Dynamic Treatment Effects (FEct Approach)

plot(fect_result, main = "Estimated ATT", 
     ylab = "Effect on Democratic Performance",
     xlab = "Years Since Treatment", 
     cex.lab = 0.9, cex.axis = 0.7, count = F,
     ylim = c(-1.5, 0.5), xlim = c(-10, 10))

## Figure 4: Diagnostic Test (FEct Approach)

plot(fect_result, type = "equiv", ylim = c(-0.25, 0.25),
     cex.legend = 0.5, cex.lab = 0.8, cex.axis = 0.7, 
     main = "Test for No Pre-Trend", stats = c("equiv.p"), 
     ylab = "Effect on Democratic Performance",
     xlab = "Years Since Treatment", count = F)

# TWFE model controlling for party

df <- df %>% 
  mutate(partycontrol_fac = as.factor(partycontrol)) %>%
  mutate(partycontrol_fac = relevel(partycontrol_fac, ref = "1"))

did_2 <- lm_robust(democracy_mcmc ~ treatment_dummy + partycontrol_fac + 
                     factor(st) + factor(year), 
                   clusters = st, data = df, se_type = "stata")

## Table 2: Effect on State Democratic Performance (Controlling for Party)

texreg(list(did_2), include.ci = FALSE, 
       omit.coef = c('st|year|Intercept'), 
       custom.coef.names = c("ATT", "Republican Party Control", "Democratic Party Control"),
       caption = "Effect on State Democratic Performance: Controlling for Party",
       custom.model.names = c("Model 2"),
       custom.gof.rows = list("State/Year Fixed Effects" = c("Yes")))

# Appendices ----

# load cost of voting data (lower values are better)

covi_data <- read.csv("covi.csv")

covi_data <- covi_data[, c('state','year', 'FinalCOVI')]

df2 <- left_join(df, covi_data, by = c("st" = "state", "year"))

# rescale to mean = 0, sd = 1

df2$scaledCOVI <- scale(df2$FinalCOVI)

# TWFE model with COVI outcome measure

did_covi <- lm_robust(scaledCOVI ~ treatment_dummy + factor(st) + factor(year),
                      clusters = st, data = df2, se_type = "stata")

# load gerrymandering data (positive values indicate Democratic advantage)

load("state_house_analysis.RData")
load("congress_analysis.RData")

state_elections_statehouse$state <- state_elections_statehouse$stateabrev
state_elections_statehouse$year <- state_elections_statehouse$cycle
state_elections_congress$state <- state_elections_congress$STATEAB
state_elections_congress$year <- state_elections_congress$election
gerrymander_st <- state_elections_statehouse[, c("year", "state", "factor")]
gerrymander_con <- state_elections_congress[, c("year", "state", "factor")]

df3 <- left_join(df, gerrymander_st, by = c("st" = "state", "year"))
df4 <- left_join(df, gerrymander_con, by = c("st" = "state", "year"))

# rescale to mean = 0, sd = 1

df3$scaled_factor_st <- scale(df3$factor)
df4$scaled_factor_con <- scale(df4$factor)

# TWFE models with gerrymandering outcome measures

did_gm_st <- lm_robust(scaled_factor_st ~ treatment_dummy + factor(st) + factor(year),
                       clusters = st, data = df3, se_type = "stata")

did_gm_con <- lm_robust(scaled_factor_con ~ treatment_dummy + factor(st) + factor(year),
                        clusters = st, data = df4, se_type = "stata")

## Table B1: Effect on State Democratic Performance (Alternative Outcome Measures)

texreg(list(did_covi, did_gm_con, did_gm_st), include.ci = FALSE, 
       omit.coef = c('st|year|Intercept'), 
       custom.coef.names = c("ATT"),
       caption = "Effect on State Democratic Performance: Alternative Outcome Measures",
       custom.model.names = c("Cost of Voting Index", "Partisan Advantage in
                              Congressional Districting", "Partisan Advantage in State
                              House Districting"),
       custom.gof.rows = list("State/Year Fixed Effects" = c("Yes", "Yes", "Yes")))

# controlled direct effects approach

set.seed(2010)

df <- df %>%
  mutate(rep_control = ifelse(partycontrol == 0, 1, 0))

direct <- sequential_g(democracy_mcmc ~ treatment_dummy + factor(st) + factor(year) | factor(st) + factor(year) | rep_control, data = df)

summary(direct)

# mediation analysis 

set.seed(2010)

med.fit <- lm(rep_control ~ treatment_dummy + factor(st) + factor(year), data = df)
out.fit <- lm(democracy_mcmc ~ rep_control + treatment_dummy + factor(st) + factor(year), data = df)

med.out <- mediate(med.fit, out.fit, treat = "treatment_dummy", 
                   mediator = "rep_control",
                   robustSE = T, sims = 1000)

summary(med.out)

## Figure C1: Causal Mediation Analysis

plot(med.out)

# freezing party control

df <- df %>%
  mutate(partycontrol_frozen = partycontrol)

for(i in levels(factor(df$state))){
  df$partycontrol_frozen[df$state == i & df$year > 2010] <- df$partycontrol[df$state == i & df$year == 2010]
}

df <- df %>% 
  mutate(partycontrol_frozen_fac = as.factor(partycontrol_frozen)) %>%
  mutate(partycontrol_frozen_fac = relevel(partycontrol_frozen_fac, ref = "1"))

# TWFE model with frozen party control

did_3 <- lm_robust(democracy_mcmc ~ treatment_dummy + partycontrol_frozen_fac +
                     factor(st) + factor(year), 
                   clusters = st, data = df, se_type = "stata")

## Table D1: Effect on State Democratic Performance (Controlling for Pre-2010 Party)

texreg(list(did_3), include.ci = FALSE, 
       omit.coef = c('st|year|Intercept'), 
       custom.coef.names = c("ATT", "Republican Party Control (Prior to CU)", "Democratic Party Control (Prior to CU)"),
       caption = "Effect on State Democratic Performance: Controlling for Party (Pre-2010)",
       custom.model.names = c("Model 3"),
       custom.gof.rows = list("State/Year Fixed Effects" = c("Yes")))


